Fit condition model with environmental and biological predictors

Fit main model (see exploratory scripts and model comparison), visualize results.

Load packages

library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(viridis)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
#> Warning: package 'sf' was built under R version 4.0.5
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(gganimate)
library(gifski)
library(latex2exp)
library(patchwork)
library(png)
library(RCurl)
library(wesanderson)
library(qwraps2) 
library(ggcorrplot)
library(ggridges)
library(brms)
#> Warning: package 'Rcpp' was built under R version 4.0.5
# remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB) # sdmTMB_0.1.4.9004

# To load entire cache in interactive r session, do: 
# qwraps2::lazyload_cache_dir(path = "R/analysis/condition_model_cf_cache/html")

Code for map plots

# Specify map ranges
ymin = 52; ymax = 60; xmin = 10; xmax = 24

map_data <- rnaturalearth::ne_countries(
  scale = "medium",
  returnclass = "sf", continent = "europe")

# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
sf::sf_use_s2(FALSE)
# https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data

swe_coast <- suppressWarnings(suppressMessages(
  st_crop(map_data,
          c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))

# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)

# Define plotting theme for main plot
theme_plot <- function(base_size = 11, base_family = "") {
  theme_light(base_size = base_size, base_family = "") +
    theme(
      legend.position = "bottom",
      legend.key.height = unit(0.2, "cm"),
      legend.margin = margin(0, 0, 0, 0),
      legend.box.margin = margin(-5, -5, -5, -5),
      strip.text = element_text(size = 9, colour = 'gray10', margin = margin(b = 1, t = 1)),
      strip.background = element_rect(fill = "grey95")
      )
}

# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 11, base_family = "") {
  theme_light(base_size = base_size, base_family = "") +
    theme(
        axis.text.x = element_text(angle = 90),
        axis.text = element_text(size = 7),
        strip.text = element_text(size = 8, colour = 'gray10', margin = margin(b = 1, t = 1)),
        strip.background = element_rect(fill = "gray95"),
        legend.direction = "horizontal",
        legend.margin = margin(1, 1, 1, 1),
        legend.box.margin = margin(0, 0, 0, 0),
        legend.key.height = unit(0.4, "line"),
        legend.key.width = unit(2, "line"),
        legend.spacing.x = unit(0.1, 'cm'),
        legend.position = "bottom",
      )
}

# Make default base map plot
xmin2 <- 303000; xmax2 <- 916000; xrange <- xmax2 - xmin2
ymin2 <- 5980000; ymax2 <- 6450000; yrange <- ymax2 - ymin2

plot_map_raster <- 
ggplot(swe_coast_proj) + 
  xlim(xmin2, xmax2) +
  ylim(ymin2, ymax2) +
  labs(x = "Longitude", y = "Latitude") +
  geom_sf(size = 0.3) + 
  theme_plot() +
  guides(colour = guide_colorbar(title.position = "top", title.hjust = 0.5),
         fill = guide_colorbar(title.position = "top", title.hjust = 0.5))

Read data

# Read the split files and join them
d1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod-condition/master/data/for_analysis/mdat_cond_(1_2).csv")
d2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod-condition/master/data/for_analysis/mdat_cond_(2_2).csv")

d <- bind_rows(d1, d2)

unique(is.na(d$density_cod))
#> [1] FALSE
unique(is.na(d$density_cod_rec))
#> [1] FALSE  TRUE

# Calculate standardized variables
d <- d %>% 
  mutate(ln_length_cm = log(length_cm),
         ln_weight_g = log(weight_g),
         year_ct = year - mean(year),
         biomass_her_sc = biomass_her,
         biomass_her_sd_sc = biomass_her_sd,
         biomass_spr_sc = biomass_spr,
         biomass_spr_sd_sc = biomass_spr_sd,
         density_cod_sc = density_cod,
         density_cod_rec_sc = density_cod_rec,
         density_fle_sc = density_fle,
         density_fle_rec_sc = density_fle_rec,
         density_saduria_sc = density_saduria,
         density_saduria_rec_sc = density_saduria_rec,
         depth_sc = depth,
         depth_rec_sc = depth_rec,
         oxy_sc = oxy,
         oxy_rec_sc = oxy_rec,
         temp_sc = temp,
         temp_rec_sc = temp_rec,
         year_f = as.factor(year))  %>%
  mutate_at(c("biomass_her_sc", "biomass_her_sd_sc", "biomass_spr_sc", "biomass_spr_sd_sc",
              "density_cod_sc", "density_cod_rec_sc", 
              "density_fle_sc", "density_fle_rec_sc", 
              "density_saduria_sc", "density_saduria_rec_sc", 
              "depth_sc", "depth_rec_sc", 
              "oxy_sc", "oxy_rec_sc", 
              "temp_sc", "temp_rec_sc"
              ),
            ~(scale(.) %>% as.vector)) %>% 
  mutate(year = as.integer(year))

unique(is.na(d))
#>      weight_g length_cm  year   lat   lon ices_rect sub_div depth   oxy  temp
#> [1,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [2,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [3,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [4,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [5,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [6,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [7,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#>          X     Y density_cod density_fle density_cod_data density_fle_data
#> [1,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [2,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [3,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [4,] FALSE FALSE       FALSE       FALSE             TRUE             TRUE
#> [5,] FALSE FALSE       FALSE       FALSE             TRUE             TRUE
#> [6,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [7,] FALSE FALSE       FALSE       FALSE             TRUE             TRUE
#>      density_saduria density_saduria_rec density_cod_rec density_fle_rec
#> [1,]           FALSE               FALSE           FALSE           FALSE
#> [2,]           FALSE               FALSE           FALSE           FALSE
#> [3,]           FALSE               FALSE           FALSE           FALSE
#> [4,]           FALSE               FALSE           FALSE           FALSE
#> [5,]           FALSE                TRUE            TRUE            TRUE
#> [6,]           FALSE                TRUE            TRUE            TRUE
#> [7,]           FALSE               FALSE           FALSE           FALSE
#>      depth_rec oxy_rec temp_rec biomass_spr biomass_her biomass_spr_sd
#> [1,]     FALSE   FALSE    FALSE       FALSE       FALSE          FALSE
#> [2,]     FALSE   FALSE    FALSE        TRUE        TRUE          FALSE
#> [3,]     FALSE   FALSE    FALSE        TRUE        TRUE           TRUE
#> [4,]     FALSE   FALSE    FALSE       FALSE       FALSE          FALSE
#> [5,]      TRUE    TRUE     TRUE        TRUE        TRUE          FALSE
#> [6,]      TRUE    TRUE     TRUE        TRUE        TRUE          FALSE
#> [7,]     FALSE   FALSE    FALSE        TRUE        TRUE          FALSE
#>      biomass_her_sd ln_length_cm ln_weight_g year_ct biomass_her_sc
#> [1,]          FALSE        FALSE       FALSE   FALSE          FALSE
#> [2,]          FALSE        FALSE       FALSE   FALSE           TRUE
#> [3,]           TRUE        FALSE       FALSE   FALSE           TRUE
#> [4,]          FALSE        FALSE       FALSE   FALSE          FALSE
#> [5,]          FALSE        FALSE       FALSE   FALSE           TRUE
#> [6,]          FALSE        FALSE       FALSE   FALSE           TRUE
#> [7,]          FALSE        FALSE       FALSE   FALSE           TRUE
#>      biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc density_cod_sc
#> [1,]             FALSE          FALSE             FALSE          FALSE
#> [2,]             FALSE           TRUE             FALSE          FALSE
#> [3,]              TRUE           TRUE              TRUE          FALSE
#> [4,]             FALSE          FALSE             FALSE          FALSE
#> [5,]             FALSE           TRUE             FALSE          FALSE
#> [6,]             FALSE           TRUE             FALSE          FALSE
#> [7,]             FALSE           TRUE             FALSE          FALSE
#>      density_cod_rec_sc density_fle_sc density_fle_rec_sc density_saduria_sc
#> [1,]              FALSE          FALSE              FALSE              FALSE
#> [2,]              FALSE          FALSE              FALSE              FALSE
#> [3,]              FALSE          FALSE              FALSE              FALSE
#> [4,]              FALSE          FALSE              FALSE              FALSE
#> [5,]               TRUE          FALSE               TRUE              FALSE
#> [6,]               TRUE          FALSE               TRUE              FALSE
#> [7,]              FALSE          FALSE              FALSE              FALSE
#>      density_saduria_rec_sc depth_sc depth_rec_sc oxy_sc oxy_rec_sc temp_sc
#> [1,]                  FALSE    FALSE        FALSE  FALSE      FALSE   FALSE
#> [2,]                  FALSE    FALSE        FALSE  FALSE      FALSE   FALSE
#> [3,]                  FALSE    FALSE        FALSE  FALSE      FALSE   FALSE
#> [4,]                  FALSE    FALSE        FALSE  FALSE      FALSE   FALSE
#> [5,]                   TRUE    FALSE         TRUE  FALSE       TRUE   FALSE
#> [6,]                   TRUE    FALSE         TRUE  FALSE       TRUE   FALSE
#> [7,]                  FALSE    FALSE        FALSE  FALSE      FALSE   FALSE
#>      temp_rec_sc year_f
#> [1,]       FALSE  FALSE
#> [2,]       FALSE  FALSE
#> [3,]       FALSE  FALSE
#> [4,]       FALSE  FALSE
#> [5,]        TRUE  FALSE
#> [6,]        TRUE  FALSE
#> [7,]       FALSE  FALSE

unique(is.na(d$density_cod_rec))
#> [1] FALSE  TRUE
unique(is.na(d$density_cod))
#> [1] FALSE

# Drop NA covariates for the correlation plot
d <- d %>% drop_na(biomass_her_sc, biomass_her_sd_sc, biomass_spr_sc, biomass_spr_sd_sc,
                   density_cod_sc, density_cod_rec_sc, density_fle_sc, density_fle_rec_sc, 
                   density_saduria_sc, density_saduria_rec_sc, depth_sc, depth_rec_sc,
                   oxy_sc, oxy_rec_sc, temp_sc, temp_rec_sc)

# Sample size
nrow(d)
#> [1] 94295

# stdev of oxygen
sd(d$oxy)
#> [1] 1.805214


# Plot correlation between variables
d_cor <- d %>% dplyr::select("biomass_her_sc", "biomass_her_sd_sc",
                             "biomass_spr_sc", "biomass_spr_sd_sc",
                             "density_cod_sc", "density_cod_rec_sc", 
                             "density_fle_sc", "density_fle_rec_sc", 
                             "density_saduria_sc", "density_saduria_rec_sc", 
                             "depth_sc", "depth_rec_sc",
                             "oxy_sc", "oxy_rec_sc",
                             "temp_sc", "temp_rec_sc")

corr <- round(cor(d_cor), 1)

ggcorrplot(corr, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 2.5) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.3))

ggsave("figures/supp/condition/correlation_variables.png", width = 6.5, height = 6.5, dpi = 600)

# Plot data
d_plot <- d %>% 
  group_by(year, X, Y) %>% 
  summarise(n = n())

plot_map_raster + 
  geom_point(data = d_plot, aes(X*1000, Y*1000, size = n, color = n), alpha = 0.5) + 
  scale_color_viridis() + 
  scale_size(range = c(.1, 3), name = "N per haul") + 
  facet_wrap(~year, ncol = 5) + 
  guides(size = "none") + 
  theme_facet_map()

ggsave("figures/supp/condition/spatial_cond_data.png", width = 6.5, height = 8.5, dpi = 600)

Read the prediction grids

pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")

pred_grid <- bind_rows(pred_grid1, pred_grid2)

unique(is.na(pred_grid$density_cod))
#> [1] FALSE
unique(is.na(pred_grid$density_cod_rec))
#> [1] FALSE

pred_grid <- pred_grid %>% mutate(year = as.integer(year),
                                  year_f = as.factor(year))

# Scale the variables with respect to data! First calculate means in data
data_means <- d %>%
  dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
                density_cod, density_cod_rec, 
                density_fle, density_fle_rec, 
                density_saduria, density_saduria_rec, 
                depth, depth_rec,
                oxy, oxy_rec, 
                temp, temp_rec) %>%
  mutate_all(~mean(.)) %>%
  distinct(.keep_all = TRUE)

data_stdev <- d %>%
  dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
                density_cod, density_cod_rec, 
                density_fle, density_fle_rec, 
                density_saduria, density_saduria_rec, 
                depth, depth_rec, 
                oxy, oxy_rec, 
                temp, temp_rec
                ) %>%
  mutate_all(~sd(.)) %>%
  distinct(.keep_all = TRUE)

# Before actually scaling, replace the ices-rectangle pelagic values that are NA with the mean across rectangles in the sub division.
# Replace the sub-division pelagic values that are NA with the mean in the year.
# Note that the sub-division values are not the sum of the rectangles (due to missing rectangles), so 
# I need to calculate a sub-division mean across rectangles within each sub-division
pred_grid <- pred_grid %>% 
  group_by(year) %>% 
  mutate(median_biomass_sprat_across_sub_div = median(biomass_spr_sd, na.rm = TRUE),
         median_biomass_herring_across_sub_div = median(biomass_her_sd, na.rm = TRUE)) %>% 
  ungroup() %>% # Replace sub_divsion NA's with the median across sub_divisions in that year
  mutate(biomass_spr_sd = ifelse(is.na(biomass_spr_sd) == TRUE, median_biomass_sprat_across_sub_div, biomass_spr_sd),
         biomass_her_sd = ifelse(is.na(biomass_her_sd) == TRUE, median_biomass_herring_across_sub_div, biomass_her_sd)) %>% 
  group_by(year, sub_div) %>% 
  mutate(median_biomass_sprat_across_rect_in_sub_div = median(biomass_spr, na.rm = TRUE),
         median_biomass_herring_across_rect_in_sub_div = median(biomass_her, na.rm = TRUE)) %>% 
  ungroup() %>% # Replace rectangle NA's with the median across rectangles in that year and sub-division
  mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, median_biomass_sprat_across_rect_in_sub_div, biomass_spr),
         biomass_her = ifelse(is.na(biomass_her) == TRUE, median_biomass_herring_across_rect_in_sub_div, biomass_her)) %>% 
  # Since I still have some NAs (some sub-divisions do not have a single rectangle in some years), I will will those rectangles 
  # with the sub-division value divided by the number of rectangles in that sub division
  group_by(year, sub_div) %>% 
  mutate(n_rect = length(unique(ices_rect))) %>% 
  ungroup() %>% 
  mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, biomass_spr_sd/n_rect, biomass_spr),
         biomass_her = ifelse(is.na(biomass_her) == TRUE, biomass_her_sd/n_rect, biomass_her))

pred_grid <- pred_grid %>%
  mutate(ln_length_cm = log(1),
         year_ct = year - mean(year),
         biomass_her_sc = (biomass_her - data_means$biomass_her) / data_stdev$biomass_her,
         biomass_her_sd_sc = (biomass_her_sd - data_means$biomass_her_sd) / data_stdev$biomass_her_sd,
         biomass_spr_sc = (biomass_spr - data_means$biomass_spr) / data_stdev$biomass_spr,
         biomass_spr_sd_sc = (biomass_spr_sd - data_means$biomass_spr_sd) / data_stdev$biomass_spr_sd,
         density_cod_sc = (density_cod - data_means$density_cod) / data_stdev$density_cod,
         density_cod_rec_sc = (density_cod_rec - data_means$density_cod_rec) / data_stdev$density_cod_rec,
         density_fle_sc = (density_fle - data_means$density_fle) / data_stdev$density_fle,
         density_fle_rec_sc = (density_fle_rec - data_means$density_fle_rec) / data_stdev$density_fle_rec,
         density_saduria_sc = (density_saduria - data_means$density_saduria) / data_stdev$density_saduria,
         density_saduria_rec_sc = (density_saduria_rec - data_means$density_saduria_rec) / data_stdev$density_saduria_rec,
         depth_sc = (depth - data_means$depth) / data_stdev$depth,
         depth_rec_sc = (depth_rec - data_means$depth_rec) / data_stdev$depth_rec,
         oxy_sc = (oxy - data_means$oxy) / data_stdev$oxy,
         oxy_rec_sc = (oxy_rec - data_means$oxy_rec) / data_stdev$oxy_rec,
         temp_sc = (temp - data_means$temp) / data_stdev$temp,
         temp_rec_sc = (temp_rec - data_means$temp_rec) / data_stdev$temp_rec,
         )

# Plot on map:
ggplot(pred_grid2, aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

pred_grid %>% 
  drop_na() %>% 
  ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

pred_grid %>% 
  drop_na(biomass_spr_sd) %>% 
  ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

unique(is.na(pred_grid))
#>          X     Y depth  year    id oxy_q3   oxy temp_q3  temp   lon   lat
#> [1,] FALSE FALSE FALSE FALSE FALSE  FALSE FALSE   FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE  FALSE FALSE   FALSE FALSE FALSE FALSE
#>      sub_div ices_rect density_saduria biomass_spr biomass_her biomass_spr_sd
#> [1,]   FALSE     FALSE           FALSE       FALSE       FALSE          FALSE
#> [2,]   FALSE     FALSE           FALSE       FALSE       FALSE          FALSE
#>      biomass_her_sd density_cod density_fle depth_rec temp_rec oxy_rec
#> [1,]          FALSE       FALSE       FALSE     FALSE    FALSE   FALSE
#> [2,]          FALSE       FALSE       FALSE     FALSE    FALSE   FALSE
#>      density_cod_rec density_fle_rec density_saduria_rec year_f
#> [1,]           FALSE           FALSE               FALSE  FALSE
#> [2,]           FALSE           FALSE               FALSE  FALSE
#>      median_biomass_sprat_across_sub_div median_biomass_herring_across_sub_div
#> [1,]                               FALSE                                 FALSE
#> [2,]                               FALSE                                 FALSE
#>      median_biomass_sprat_across_rect_in_sub_div
#> [1,]                                       FALSE
#> [2,]                                        TRUE
#>      median_biomass_herring_across_rect_in_sub_div n_rect ln_length_cm year_ct
#> [1,]                                         FALSE  FALSE        FALSE   FALSE
#> [2,]                                          TRUE  FALSE        FALSE   FALSE
#>      biomass_her_sc biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc
#> [1,]          FALSE             FALSE          FALSE             FALSE
#> [2,]          FALSE             FALSE          FALSE             FALSE
#>      density_cod_sc density_cod_rec_sc density_fle_sc density_fle_rec_sc
#> [1,]          FALSE              FALSE          FALSE              FALSE
#> [2,]          FALSE              FALSE          FALSE              FALSE
#>      density_saduria_sc density_saduria_rec_sc depth_sc depth_rec_sc oxy_sc
#> [1,]              FALSE                  FALSE    FALSE        FALSE  FALSE
#> [2,]              FALSE                  FALSE    FALSE        FALSE  FALSE
#>      oxy_rec_sc temp_sc temp_rec_sc
#> [1,]      FALSE   FALSE       FALSE
#> [2,]      FALSE   FALSE       FALSE

pred_grid %>% 
  drop_na(biomass_spr_sc, biomass_spr_sd_sc, biomass_her_sc, biomass_her_sd_sc) %>% 
  ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

Make spde mesh

spde <- make_mesh(d, xy_cols = c("X", "Y"),
                  n_knots = 100,  
                  type = "kmeans", seed = 42)

# Plot and save spde
png(file = "figures/supp/condition/spde.png", units = "in", width = 6.5, height = 6.5, res = 300)
plot(spde)
dev.off()

Fit models

Non spatial model to get the average length weight relationship. The main model’s response variable is the ratio of the observed vs this predicted average weight.

ptm <- proc.time()
mnull <- sdmTMB(
  formula = ln_weight_g ~ ln_length_cm,
  data = d,
  time = NULL,
  mesh = spde,
  family = student(link = "identity", df = 5),
  spatiotemporal = "off",
  spatial = "off",
  silent = TRUE,
  reml = FALSE
  )
proc.time() - ptm
#>    user  system elapsed 
#>   1.611   0.116   1.763

sanity(mnull)
#> ✖ Non-linear minimizer did not converge: do not trust this model!
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# mlm <- lm(
#   formula = ln_weight_g ~ ln_length_cm,
#   data = d)

# bm <- brms::brm(formula = ln_weight_g ~ ln_length_cm, data = d, family = student())
# summary(bm)
# plot(bm)
# brms::pp_check(bm)

# Right, so even though sdmTMB warns, the actual estimates are very similar to lm and brms, student-t or gaussian. So will just play along with it

# Extract average a and b across the total dataset
a <- as.numeric(tidy(mnull)[1, "estimate"])
b <- as.numeric(tidy(mnull)[2, "estimate"])
a
#> [1] -4.638609
b
#> [1] 2.984044

d$weight_g_avg_pred <- exp(a)*d$length_cm^b

plot(d$weight_g ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$weight_g); abline(a = 0, b = 1, col = "red")

d$le_cren <- d$weight_g / d$weight_g_avg_pred
d$log_le_cren <- log(d$le_cren)
hist(d$le_cren)

Fit main model

ptm <- proc.time()
mfull <- sdmTMB(
  formula = log_le_cren ~ -1 + year_f + biomass_her_sc + biomass_her_sd_sc +
    biomass_spr_sc + biomass_spr_sd_sc +
    density_cod_sc + density_cod_rec_sc + 
    density_fle_sc + density_fle_rec_sc + 
    density_saduria_sc + density_saduria_rec_sc + 
    depth_sc + depth_rec_sc + 
    oxy_sc + oxy_rec_sc + 
    temp_sc + temp_rec_sc,
    data = d,
  time = "year",
  mesh = spde, 
  family = student(link = "identity", df = 5),
  spatiotemporal = "iid",
  spatial = "on",
  silent = TRUE,
  reml = FALSE,
  control = sdmTMBcontrol(newton_loops = 1)#,
  # priors = sdmTMBpriors(b = normal(c(rep(NA, 27), rep(0, 16)),
  #                                  c(rep(NA, 27), rep(1, 16))))
  )
system("say Just finished!")
proc.time() - ptm
#>    user  system elapsed 
#> 883.707  43.842 944.490

sanity(mfull)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large

tidy(mfull, conf.int = TRUE) %>% arrange(desc(abs(estimate)))
tidy(mfull, effects = "ran_par", conf.int = TRUE) 
# Arrange coefficients
tidy(mfull, conf.int = TRUE) %>% filter(!grepl("year", term)) %>% arrange(desc(estimate))
#> filter: removed 27 rows (63%), 16 rows remaining

tidy(mfull, conf.int = TRUE) %>% filter(!grepl("year", term)) %>% arrange(estimate)
#> filter: removed 27 rows (63%), 16 rows remaining

# Average magnitude of fixed effects?
tidy(mfull, conf.int = TRUE) %>%
  filter(!grepl("year", term)) %>%
  mutate(abs_coef = abs(estimate)) %>% 
  summarise(mean_abs_coef = mean(abs_coef))
#> filter: removed 27 rows (63%), 16 rows remaining
#> mutate: new variable 'abs_coef' (double) with 16 unique values and 0% NA
#> summarise: now one row and one column, ungrouped

tidy(mfull, conf.int = TRUE, effects = "ran_par") %>%
  filter(term %in% c("sigma_O", "sigma_E")) %>%
  summarise(mean_estimate = mean(estimate))
#> filter: removed 2 rows (50%), 2 rows remaining
#> summarise: now one row and one column, ungrouped
# This model is for the sensitivity analysis, comparing estimates between different model subsets etc. These tend to fit better with only a spatiotemporal term (spatial stdev tends to blow up), so I want them to be as comparable as possible
ptm <- proc.time()
mfull_spatiotemporal <- sdmTMB(
  formula = log_le_cren ~ -1 + year_f + biomass_her_sc + biomass_her_sd_sc +
    biomass_spr_sc + biomass_spr_sd_sc +
    density_cod_sc + density_cod_rec_sc + 
    density_fle_sc + density_fle_rec_sc + 
    density_saduria_sc + density_saduria_rec_sc + 
    depth_sc + depth_rec_sc + 
    oxy_sc + oxy_rec_sc + 
    temp_sc + temp_rec_sc,
    data = d,
  time = "year",
  mesh = spde, 
  family = student(link = "identity", df = 5),
  spatiotemporal = "iid",
  spatial = "off",
  silent = TRUE,
  reml = FALSE,
  control = sdmTMBcontrol(newton_loops = 1)
  )
system("say Just finished!")
proc.time() - ptm
#>    user  system elapsed 
#> 672.519  20.072 704.688

sanity(mfull_spatiotemporal)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large

tidy(mfull_spatiotemporal, conf.int = TRUE) %>% arrange(desc(abs(estimate)))
tidy(mfull_spatiotemporal, effects = "ran_par", conf.int = TRUE) 


# Extract random and fixed coefficients from the full model
mfull_est_st <- bind_rows(tidy(mfull_spatiotemporal, effects = "ran_par", conf.int = TRUE) %>%
                            filter(term == "sigma_E"),
                          tidy(mfull_spatiotemporal, effects = "fixed", conf.int = TRUE)  %>% 
                            filter(!grepl('year', term))) %>%
  
  mutate(term = factor(term))
#> filter: removed 2 rows (67%), one row remaining
#> filter: removed 27 rows (63%), 16 rows remaining
#> mutate: converted 'term' from character to factor (0 new NA)

# Extract coefficients and save as csv
mfull_est_st <- mfull_est_st %>% 
  mutate(group = "Herring",
         group = ifelse(term %in% c("biomass_spr_sc", "biomass_spr_sd_sc"), "Sprat", group),
         group = ifelse(term %in% c("density_cod_rec_sc", "density_cod_sc"), "Cod", group),
         group = ifelse(term %in% c("density_fle_rec_sc", "density_fle_sc"), "Flounder", group),
         group = ifelse(term %in% c("density_saduria_rec_sc", "density_saduria_sc"), "Saduria", group),
         group = ifelse(term %in% c("depth_rec_sc", "depth_sc"), "Depth", group),
         group = ifelse(term %in% c("oxy_rec_sc", "oxy_sc"), "Oxygen", group),
         group = ifelse(term %in% c("temp_rec_sc", "temp_sc"), "Temp", group),
         group = ifelse(term == "sigma_E", "Random", group),
         ) %>% 
  mutate(scale = "Subdivision",
         scale = ifelse(term %in% c("biomass_her_sc", "biomass_spr_sc", "density_cod_rec_sc", "density_fle_rec_sc",
                                    "density_saduria_rec_sc", "depth_rec_sc", "oxy_rec_sc", "temp_rec_sc"), "Rectangle", scale),
         scale = ifelse(term %in% c("density_cod_sc", "density_fle_sc", "density_saduria_sc",
                                    "depth_sc", "oxy_sc", "temp_sc"), "Haul", scale),
         scale = ifelse(term == "sigma_E", "Spatial/spatiotemporal s.d.", scale),
         ) %>% 
  mutate(term2 = ifelse(term == "sigma_E", 2, 1))
#> mutate: new variable 'group' (character) with 9 unique values and 0% NA
#> mutate: new variable 'scale' (character) with 4 unique values and 0% NA
#> mutate: new variable 'term2' (double) with 2 unique values and 0% NA

# Now save this for comparison in sensitivity analysis script
write.csv(mfull_est_st, "output/mfull_est_st.csv")
# Save model (so that I can add density predictions to the pred grid for condition)
saveRDS(mfull, "output/mfull.rds")
mcmc_res <- residuals(mfull, type = "mle-mcmc", mcmc_iter = 5000, mcmc_warmup = 2000)

png(file = "figures/supp/condition/qq.png", units = "in", width = 6.5, height = 6.5, res = 300)
qqnorm(mcmc_res);qqline(mcmc_res)
dev.off()

d$residuals <- mcmc_res[, 1]

Calculate approximate \(R^2\) manually using a modifed r2.sdmTMB

b_fe <- tidy(mfull)
b_fe <- stats::setNames(b_fe$estimate, b_fe$term)
X <- mfull$tmb_data$X_ij
X <- matrix(unlist(X), ncol = length(b_fe))

VarF <- var(as.vector(b_fe %*% t(X))) # variance from fixed-effects
 
# What if I skip the year effect?
b_fe2 <- tidy(mfull) %>% filter(!grepl('year', term))
#> filter: removed 27 rows (63%), 16 rows remaining
b_fe2 <- stats::setNames(b_fe2$estimate, b_fe2$term)
X2 <- mfull$tmb_data$X_ij

first_index <- match('year_f2019', names(as.data.frame(mfull$tmb_data$X_ij[[1]]))) + 1
last_index <- ncol(mfull$tmb_data$X_ij[[1]])

X2 <- matrix(unlist(X2), ncol = length(b_fe))[, first_index:last_index] # skip year columns
VarF2 <- var(as.vector(b_fe2 %*% t(X2))) # variance from fixed-effects

b_ran <- tidy(mfull, "ran_par")
sigma_O <- b_ran$estimate[b_ran$term == "sigma_O"] # spatial standard deviation
sigma_E <- b_ran$estimate[b_ran$term == "sigma_E"] # spatiotemporal standard deviation
VarO <- sigma_O^2 # spatial variance
VarE <- sigma_E^2 # spatiotemporal variance
 
# Calculate obs - pred
d_r2 <- d
d_r2$pred <- predict(mfull)
d_r2$resid_manual <- d_r2$le_cren - d_r2$pred$est
#hist(d_r2$resid_manual)

VarR <- var(as.vector(d_r2$resid_manual)) # "residual" variance

# Marginal R2s
denominator <- VarF + VarO + VarE + VarR
denominator_no_yr <- VarF2 + VarO + VarE + VarR

marg <- VarF/denominator
marg_no_yr <- VarF2/denominator_no_yr

out <- data.frame(
    marginal = marg,
    partial_spatial = VarO/denominator,
    partial_spatiotemporal = VarE/denominator,
    conditional = (VarF + VarO + VarE)/denominator,
    all_random = (VarO + VarE)/denominator,
    marginal_random_ratio = marg / ((VarO + VarE)/denominator),
    random_marginal_ratio = ((VarO + VarE)/denominator) / marg
  )
out

out_no_yr <- data.frame(
    marginal = marg_no_yr,
    partial_spatial = VarO/denominator_no_yr,
    partial_spatiotemporal = VarE/denominator_no_yr,
    conditional = (VarF2 + VarO + VarE)/denominator_no_yr,
    all_random = (VarO + VarE)/denominator_no_yr,
    marginal_random_ratio = marg_no_yr / ((VarO + VarE)/denominator_no_yr),
    random_marginal_ratio = ((VarO + VarE)/denominator_no_yr) / marg_no_yr
  )
out_no_yr

Plot residuals

# Residuals on map
plot_map_raster +
  geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = residuals), size = 0.5) +
  scale_colour_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  labs(color = "residuals") +
  theme_facet_map() +
  geom_sf(size = 0.1)


ggsave("figures/supp/condition/spatial_resid.png", width = 6.5, height = 8.5, dpi = 600)

# Residuals vs length and year
ggplot(d, aes(x = length_cm, y = residuals, color = residuals)) +
  geom_point(size = 0.1, alpha = 1) +
  geom_hline(yintercept = 0, size = 0.5, color = "black", linetype = 2, alpha = 0.5) +
  labs(x = "length [cm]", y = "Residuals", color = "") +
  stat_smooth(color = "black", size = 0.5) +
  facet_wrap(~year, ncol = 5) +
  scale_color_gradient2() +
  theme_plot()


ggsave("figures/supp/condition/residuals_vs_length_and_year.png", width = 6.5, height = 8.5, dpi = 600)

ggplot(d, aes(year, length_cm, z = residuals)) +
  stat_summary_2d(bins = 40) +
  scale_fill_gradient2()


# ggsave("figures/supp/condition/residuals_vs_length_and_year2.png", width = 6.5, height = 6.5, dpi = 600)

Annual index using get_index_sims

pred_avg_sim <- predict(mfull, newdata = pred_grid %>% filter(depth > 10 & depth < 110), nsim = 500)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining

#ncells <- filter(pred_grid, year == max(pred_grid$year)) %>% nrow()
# index_avg_sim <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)))
# index_avg_sims <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)),
#                                  return_sims = TRUE) # This is for calculation on samples!

# Old version above: now we do a weighted version. We could do it directly in the get_index_sims or manually in the prediction grid (but then we don't get an CIs)
# This is so that we can divide by the sum of weights (needs to be done by year)

weight_sum <- pred_grid %>% 
  filter(depth > 10 & depth < 110) %>% 
  group_by(year) %>% 
  summarise(density_cod = sum(density_cod))
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped

index_avg_sim <- get_index_sims(pred_avg_sim,
                                area = filter(pred_grid, depth > 10 & depth < 110)$density_cod) %>% 
  left_join(weight_sum) %>% 
  mutate(
  est = est / density_cod,
  lwr = lwr / density_cod,
  upr = upr / density_cod
  )
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining, by = "year"
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 27
#> > ====
#> > rows total 27
#> mutate: changed 27 values (100%) of 'est' (0 new NA)
#> changed 27 values (100%) of 'lwr' (0 new NA)
#> changed 27 values (100%) of 'upr' (0 new NA)

index_avg_sims <- get_index_sims(pred_avg_sim,
                                 area = filter(pred_grid, depth > 10 & depth < 110)$density_cod,
                                 return_sims = TRUE) %>% 
  left_join(weight_sum) %>% 
  mutate(est = .value / density_cod)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining, by = "year"
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 13,500
#> > ========
#> > rows total 13,500
#> mutate: new variable 'est' (double) with 13,500 unique values and 0% NA

ggplot(index_avg_sim, aes(year, est, ymin = lwr, ymax = upr)) +
  geom_ribbon(alpha = 0.2, color = NA) +
  geom_line(size = 1) + 
  geom_line(data = filter(index_avg_sims, .iteration < 26),
            aes(year, est, group = .iteration), inherit.aes = FALSE, alpha = 0.3)
#> filter: removed 12,825 rows (95%), 675 rows remaining

Get sims from prediction with only fixed effects (-year) and random effects

pred_grid_fe <- pred_grid %>% 
  mutate(year_f = as.factor(1993))
#> mutate: changed 249,002 values (96%) of 'year_f' (0 new NA)
  
pred_avg_sim_fe <- predict(mfull, newdata = pred_grid_fe %>% filter(depth > 10 & depth < 110), nsim = 500, re_form = NA)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining

weight_sum <- pred_grid %>% 
  filter(depth > 10 & depth < 110) %>% 
  group_by(year) %>% 
  summarise(density_cod = sum(density_cod))
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped

index_avg_sim_fe <- get_index_sims(pred_avg_sim_fe,
                                   area = filter(pred_grid, depth > 10 & depth < 110)$density_cod) %>% 
  left_join(weight_sum) %>% 
  mutate(
  est = est / density_cod,
  lwr = lwr / density_cod,
  upr = upr / density_cod
  )
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining, by = "year"
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 27
#> > ====
#> > rows total 27
#> mutate: changed 27 values (100%) of 'est' (0 new NA)
#> changed 27 values (100%) of 'lwr' (0 new NA)
#> changed 27 values (100%) of 'upr' (0 new NA)

index_avg_sims_fe <- get_index_sims(pred_avg_sim_fe,
                                    area = filter(pred_grid, depth > 10 & depth < 110)$density_cod,
                                    return_sims = TRUE) %>% 
  left_join(weight_sum) %>% 
  mutate(est = .value / density_cod)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining, by = "year"
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 13,500
#> > ========
#> > rows total 13,500
#> mutate: new variable 'est' (double) with 13,500 unique values and 0% NA

# With random effects but no fixed year effect
# pred_grid_re <- pred_grid %>% 
#   mutate(year_f = as.factor(1993))
#   
# pred_avg_sim_re <- predict(mfull, newdata = pred_grid_re %>% filter(depth > 10 & depth < 110), nsim = 500)
# 
# weight_sum <- pred_grid %>% 
#   filter(depth > 10 & depth < 110) %>% 
#   group_by(year) %>% 
#   summarise(density_cod = sum(density_cod))
# 
# index_avg_sim_re <- get_index_sims(pred_avg_sim_re,
#                                    area = filter(pred_grid, depth > 10 & depth < 110)$density_cod) %>% 
#   left_join(weight_sum) %>% 
#   mutate(
#   est = est / density_cod,
#   lwr = lwr / density_cod,
#   upr = upr / density_cod
#   )
# 
# index_avg_sims_re <- get_index_sims(pred_avg_sim_re,
#                                     area = filter(pred_grid, depth > 10 & depth < 110)$density_cod,
#                                     return_sims = TRUE) %>% 
#   left_join(weight_sum) %>% 
#   mutate(est = .value / density_cod)

# Plot all together
index_comp_sim <- bind_rows(#index_avg_sim_re %>% mutate(Prediction = "Random & Fixed effects (no yr)"),
                            index_avg_sim_fe %>% mutate(Prediction = "Fixed effects (no yr)"),
                            index_avg_sim %>% mutate(Prediction = "Full"))
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA

index_comp_sims <- bind_rows(#index_avg_sims_re %>% mutate(Prediction = "Random & Fixed effects (no yr)"),
                             index_avg_sims_fe %>% mutate(Prediction = "Fixed effects (no yr)"),
                             index_avg_sims %>% mutate(Prediction = "Full")) %>%
  filter(.iteration < 26)
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
#> filter: removed 25,650 rows (95%), 1,350 rows remaining

ggplot(index_comp_sim, aes(year, est, ymin = lwr, ymax = upr, color = Prediction, fill = Prediction)) +
  geom_ribbon(alpha = 0.2, color = NA) +
  geom_line(size = 1) + 
  geom_line(data = filter(index_comp_sims, .iteration < 26),
            aes(year, est, group = interaction(.iteration, Prediction), color = Prediction), inherit.aes = FALSE, alpha = 0.3) +
  NULL
#> filter: no rows removed

Difference in Le Cren condition factor between the end and start of time period using sims

# Reshape data again to calculate row-wise differences
year_diff <- index_avg_sims %>% 
  dplyr::select(-.value, -density_cod) %>% 
  filter(year %in% c(1993, 2019)) %>% 
  pivot_wider(names_from = year, values_from = est) %>% 
  mutate("2019-1993" = `2019` - `1993`) %>% 
  pivot_longer(2:4) %>% 
  ungroup()
#> filter: removed 12,500 rows (93%), 1,000 rows remaining
#> pivot_wider: reorganized (year, est) into (1993, 2019) [was 1000x3, now 500x3]
#> mutate: new variable '2019-1993' (double) with 500 unique values and 0% NA
#> pivot_longer: reorganized (1993, 2019, 2019-1993) into (name, value) [was 500x4, now 1500x3]
#> ungroup: no grouping variables

# Plot these
p1 <- ggplot(year_diff) +
  geom_density(data = filter(year_diff, !name == "2019-1993"),
               aes(value, fill = name), alpha = 0.6) +
  scale_fill_brewer(palette = "Dark2") +
  coord_cartesian(expand = 0)
#> filter: removed 500 rows (33%), 1,000 rows remaining

p2 <- ggplot(year_diff) +
  geom_density(data = filter(year_diff, name == "2019-1993"),
               aes(value), alpha = 0.6, fill = "grey") +
  geom_vline(xintercept = 0, linetype = 2) 
#> filter: removed 1,000 rows (67%), 500 rows remaining

p1 / p2


# Calculate quantiles for each level (to go in the results section)
year_diff %>% 
  group_by(name) %>% 
  summarise(median = quantile(value, probs = 0.5),
            upr97.5 = quantile(value, probs = 0.975),
            lwr2.5 = quantile(value, probs = 0.025))
#> group_by: one grouping variable (name)
#> summarise: now 3 rows and 4 columns, ungrouped

# Percent change in condition factor
year_diff %>% 
  group_by(name) %>% 
  pivot_wider(names_from = name, values_from = value) %>%
  mutate(perc_change = ((`2019`-`1993`)/`1993`)*100) %>% 
  ungroup() %>% 
  summarise(lwr2.5 = quantile(perc_change, probs = 0.025),
            median = quantile(perc_change, probs = 0.5),
            upr97.5 = quantile(perc_change, probs = 0.975))
#> group_by: one grouping variable (name)
#> pivot_wider: reorganized (name, value) into (1993, 2019, 2019-1993) [was 1500x3, now 500x4]
#> mutate: new variable 'perc_change' (double) with 500 unique values and 0% NA
#> ungroup: no grouping variables
#> summarise: now one row and 3 columns, ungrouped

Evaluate sensitivity to using a pred grid with a subset excluding low density areas

# Now do the same but excluding the areas not sampled in the pred grid
# ncells_sub <- filter(pred_grid, year == max(pred_grid$year) & depth > 40 & depth < 120) %>% nrow()
# pred_grid_sub <- filter(pred_grid, depth > 40 & depth < 120)
# 
# pred_avg_sim_sub <- predict(mfull, newdata = pred_grid_sub, nsim = 100)
# 
# index_avg_sim_sub <- get_index_sims(pred_avg_sim_sub, area = rep(1/ncells_sub, nrow(pred_avg_sim_sub)))
# 
# # Combine and plot & compare
# index_avg_sim_comp <- bind_rows(mutate(index_avg_sim, area = "full"),
#                                 mutate(index_avg_sim_sub, area = "subset"))
# 
# ggplot(index_avg_sim_comp, aes(year, est, ymin = lwr, ymax = upr, color = area, fill = area)) +
#   geom_line() +
#   geom_ribbon(alpha = 0.4, color = NA)

Simulated annual condition factor by sub-division

div_index_list <- list()

for(i in unique(pred_grid$sub_div)){
  
  pred_grid_sub <- pred_grid %>% filter(sub_div == i & depth > 10 & depth < 110)
  
  pred_avg_sim_sub <- predict(mfull, newdata = pred_grid_sub, nsim = 500)
  
  weight_sum_sub <- pred_grid_sub %>% 
    group_by(year) %>%
    summarise(density_cod = sum(density_cod))
  
  index_avg_sim <- get_index_sims(pred_avg_sim_sub,
                                  area = pred_grid_sub$density_cod) %>% 
    left_join(weight_sum_sub) %>% 
    mutate(est = est / density_cod,
           lwr = lwr / density_cod,
           upr = upr / density_cod) %>% 
    mutate(sub_div = i)
           
  div_index_list[[i]] <- index_avg_sim
  
}
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.

div_index <- bind_rows(div_index_list)

Condition factor over time

# Spatiotemporally standardized indicies
div_index2 <- bind_rows(div_index %>% dplyr::select(year, est, lwr, upr, sub_div) %>% mutate(sub_div = as.factor(sub_div)),
                        mutate(index_avg_sim, sub_div = as.factor("Full model prediction")))

# Add in raw mean of data:
div_index2 <- bind_rows(div_index2,
                        d %>%
                          group_by(year) %>%
                          summarise(est = exp(mean(log_le_cren))) %>%
                          mutate(sub_div = "Empirical mean"))

# Add fixed effect predictions
div_index2 <- bind_rows(div_index2,
                        index_avg_sim_fe %>% mutate(sub_div = "Fixed effects (year omitted)"))


# Add sims for total and fixed effect prediction
index_comp_sims <- bind_rows(index_avg_sims_fe %>% mutate(sub_div = "Fixed effects (year omitted)"),
                             index_avg_sims %>% mutate(sub_div = "Full model prediction")) %>%
  filter(.iteration < 26)

pal <- brewer.pal(n = 4, name = "Set1")[2:4]

# Plot!
p1 <- div_index2 %>% 
  filter(sub_div %in% c("Full model prediction", "Empirical mean", "Fixed effects (year omitted)")) %>% 
  mutate(plot_facet = "Total and data") %>% 
  mutate(plot_group = ifelse(sub_div == "Empirical mean", 1, 0)) %>% 
  ggplot(aes(year, est, color = reorder(sub_div, plot_group))) +
  labs(y = "Le Cren's condition factor", x = "", color = "") +
  geom_ribbon(aes(x = year, y = est, ymin = lwr, ymax = upr, fill = sub_div), color = NA, alpha = 0.25) +
  geom_line(data = filter(index_comp_sims, .iteration < 26),
            aes(year, est, group = interaction(.iteration, sub_div), color = sub_div), inherit.aes = FALSE, alpha = 0.3) +
  geom_line(size = 0.8) +
  facet_wrap(~plot_facet, ncol = 1) +
  scale_color_brewer(palette = "Set2") +
  scale_fill_brewer(palette = "Set2") +
  scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
  # scale_fill_manual(values = pal) +
  # scale_color_manual(values = pal) +
  theme_plot(base_size = 10) +
  theme(plot.margin = unit(c(0.4, 0.4, 0, 0.4), "cm"), 
        legend.position = c(0.22, 0.1),
        legend.direction = "horizontal",
        legend.spacing.y = unit(0, 'cm'),
        legend.key.size = unit(1, "cm"),
        legend.key.width = unit(0.5, "cm"),
        legend.title = element_text(size = 8),
        legend.text = element_text(size = 7), 
        legend.background = element_rect(fill = NA),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  guides(fill = "none",
         color = guide_legend(ncol = 1, title.position = "top", title.hjust = 0.5))

p1
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf


p2 <- div_index2 %>% 
  filter(!sub_div %in% c("Full model prediction", "Empirical mean", "Fixed effects (year omitted)")) %>% 
  mutate(plot_facet = "Subdivision") %>% 
  ggplot(aes(year, est, color = sub_div, fill = sub_div)) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  labs(y = "Le Cren's condition factor", x = "Year") +
  #geom_ribbon(aes(x = year, y = est, ymin = lwr, ymax = upr), color = NA, alpha = 0.15) + 
  geom_line() +
  facet_wrap(~plot_facet, ncol = 1) +
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2") +
  scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
  theme_plot(base_size = 10) +
  theme(plot.margin = unit(c(0, 0.4, 0.4, 0.4), "cm"), 
        legend.position = c(0.92, 0.88),
        legend.direction = "horizontal",
        legend.spacing.y = unit(0, 'cm'),
        legend.key.size = unit(0.8, "cm"),
        legend.key.width = unit(0.3, "cm"),
        legend.title = element_text(size = 8),
        legend.text = element_text(size = 7), 
        legend.background = element_rect(fill = NA),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  guides(color = guide_legend(ncol = 1, title.position = "top", title.hjust = 0.5),
         fill = "none")

p2


p1 / p2 + plot_annotation(tag_levels = "A")
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf


#ggsave("figures/cond_index.pdf", width = 11, height = 17, units = "cm")

Coefficients

# What is one standard deviation of oxygen?
sd(d$oxy)
#> [1] 1.805214

# Extract random and fixed coefficients from the full model
mfull_est <- bind_rows(tidy(mfull, effects = "ran_par", conf.int = TRUE) %>%
                         filter(term %in% c("sigma_O", "sigma_E")),
                       
                       tidy(mfull, effects = "fixed", conf.int = TRUE)  %>% 
                         filter(!grepl('year', term))) %>%
  
  mutate(term = factor(term))

mfull_est <- mfull_est %>% 
  mutate(group = "Herring",
         group = ifelse(term %in% c("biomass_spr_sc", "biomass_spr_sd_sc"), "Sprat", group),
         group = ifelse(term %in% c("density_cod_rec_sc", "density_cod_sc"), "Cod", group),
         group = ifelse(term %in% c("density_fle_rec_sc", "density_fle_sc"), "Flounder", group),
         group = ifelse(term %in% c("density_saduria_rec_sc", "density_saduria_sc"), "Saduria", group),
         group = ifelse(term %in% c("depth_rec_sc", "depth_sc"), "Depth", group),
         group = ifelse(term %in% c("oxy_rec_sc", "oxy_sc"), "Oxygen", group),
         group = ifelse(term %in% c("temp_rec_sc", "temp_sc"), "Temp", group),
         group = ifelse(term == "sigma_E", "Random", group),
         group = ifelse(term == "sigma_O", "Random", group)
         )

mfull_est <- mfull_est %>% 
  mutate(scale = "Subdivision",
         scale = ifelse(term %in% c("biomass_her_sc", "biomass_spr_sc", "density_cod_rec_sc", "density_fle_rec_sc",
                                    "density_saduria_rec_sc", "depth_rec_sc", "oxy_rec_sc", "temp_rec_sc"), "Rectangle", scale),
         scale = ifelse(term %in% c("density_cod_sc", "density_fle_sc", "density_saduria_sc",
                                    "depth_sc", "oxy_sc", "temp_sc"), "Haul", scale),
         scale = ifelse(term == "sigma_E", "Spatial/spatiotemporal s.d.", scale),
         scale = ifelse(term == "sigma_O", "Spatial/spatiotemporal s.d.", scale)
         )

# Quick calculation:
mfull_est %>% 
  mutate(par_type = ifelse(term %in% c("sigma_O", "sigma_E"), "re", "fe")) %>% 
  group_by(par_type) %>% 
  summarise(mean_est = mean(estimate)) %>% 
  pivot_wider(names_from = par_type, values_from = mean_est) %>% 
  mutate(ratio = re/fe)

# Plot effects
# This is the order changed labels should go in!
levels(mfull_est$term)
#>  [1] "biomass_her_sc"         "biomass_her_sd_sc"      "biomass_spr_sc"        
#>  [4] "biomass_spr_sd_sc"      "density_cod_rec_sc"     "density_cod_sc"        
#>  [7] "density_fle_rec_sc"     "density_fle_sc"         "density_saduria_rec_sc"
#> [10] "density_saduria_sc"     "depth_rec_sc"           "depth_sc"              
#> [13] "oxy_rec_sc"             "oxy_sc"                 "sigma_E"               
#> [16] "sigma_O"                "temp_rec_sc"            "temp_sc"

# I want the random effects to be dark gray...
pal <- brewer.pal(n = 8, name = "Dark2")
pal2 <- c(pal[1:5], "black", pal[6:8])

# Sort the terms so that random effects are at the top...
mfull_est <- mfull_est %>% 
  mutate(term2 = ifelse(term %in% c("sigma_E", "sigma_O"), 2, 1))

ggplot(mfull_est, aes(reorder(term, term2), estimate, color = group, fill = group, shape = scale)) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray40", alpha = 0.5) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
  geom_point(size = 2.5) +
  scale_color_manual(values = pal2, name = "") +
  scale_shape_manual(values = c(15, 19, 21, 17), name = "") +
  scale_fill_manual(values = rep("white", 9), name = "") +
  labs(y = "Estimate", x = "Standardized coefficient") +
  scale_x_discrete(breaks = levels(mfull_est$term),
                   labels = c(expression(Herring[rec]),
                              expression(Herring[sub]),
                              expression(Sprat[rec]),
                              expression(Sprat[sub]),
                              expression(Cod[rec]),
                              expression(Cod[haul]),
                              expression(Flounder[rec]),
                              expression(Flounder[haul]),
                              expression(Saduria~entomon[rec]),
                              expression(Saduria~entomon[haul]),
                              expression(Depth[rec]),
                              expression(Depth[haul]),
                              expression(Oxygen[rec]),
                              expression(Oxygen[haul]),
                              expression(sigma[E]),
                              expression(sigma[O]),
                              expression(Temp[rec]),
                              expression(Temp[haul])
                              )) +
  coord_flip() +
  theme_plot() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.margin = unit(c(0.2, 0.2, 0.4, 0.2), "cm")) + 
  guides(color = "none", fill = "none", shape = guide_legend(ncol = 4))

  
ggsave("figures/effect_size.pdf", width = 17, height = 17, units = "cm")

# Just a test to see the labels were alright
ggplot(mfull_est, aes(reorder(term, estimate), estimate)) +
ggplot(mfull_est, aes(term, estimate)) +
  geom_hline(yintercept = 0, linetype = 2, color = "red", alpha = 0.5) +
  geom_point(size = 2) +
  labs(x = "", y = "Standardized coefficient") +
  coord_flip()


# Now save this for comparison in sensitivity analysis script
write.csv(mfull_est, "output/mfull_est.csv")

Predictions in space

pred_map <- predict(mfull, newdata = pred_grid)

pal <- wes_palette("Zissou1", 21, type = "continuous")

plot_map_raster +
  geom_raster(data = filter(pred_map, depth < 120 & depth > 10), aes(x = X*1000, y = Y*1000, fill = exp(est))) +
  scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") +
  facet_wrap(~ year, ncol = 5) +
  theme_facet_map() +
  geom_sf(size = 0.1)


ggsave("figures/supp/condition/all_years_condition_map_covars.png", width = 6.5, height = 8.5, dpi = 600)

# And a smaller map for selected years
lab_size <- 1.7

plot_map_raster +
  geom_raster(data = filter(pred_map, year %in% c(1994, 2001, 2008, 2018) & depth < 120 & depth > 10),
              aes(x = X*1000, y = Y*1000, fill = exp(est))) +
  scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") + 
  facet_wrap(~ year, ncol = 2) +
  theme(legend.key.height = unit(0.1, "cm"),
        legend.key.width = unit(0.5, "cm"),
        legend.position = c(0.9, .036),
        legend.direction = "horizontal",
        legend.background = element_rect(fill = NA),
        legend.text = element_text(size = 6),
        legend.title = element_text(size = 7.4)) +
  geom_sf(size = 0.1) + 
  annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)


ggsave("figures/condition_map_covars.pdf", width = 17, height = 17, units = "cm")

Plot spatial and spatiotemporal random effects:

# Plot spatiotemporal random effect
plot_map_raster +
  geom_raster(data = pred_map, aes(x = X * 1000, y = Y * 1000, fill = epsilon_st)) +
  scale_fill_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  theme_facet_map() +
  geom_sf(size = 0.1)


ggsave("figures/supp/condition/epsilon_st_map.png", width = 6.5, height = 8.5, dpi = 600)

# Plot spatial random effect
plot_map_raster +
  geom_raster(data = filter(pred_map, year == 1999), aes(x = X * 1000, y = Y * 1000, fill = omega_s)) +
  scale_fill_gradient2() +
  geom_sf(size = 0.1)


ggsave("figures/supp/condition/omega_s_map.png", width = 6.5, height = 6.5, dpi = 600)

Calculate the “spatial trend” from the estimates:

# Fit a linear model to each prediction grid of the estimate over time
# https://community.rstudio.com/t/extract-slopes-by-group-broom-dplyr/2751/7
time_slopes_by_year <- pred_map %>%
  mutate(id = paste(X, Y, sep = "_")) %>%
  split(.$id) %>%
  purrr::map(~lm(est ~ year, data = .x)) %>%
  purrr::map_df(broom::tidy, .id = 'id') %>%
  filter(term == 'year')

# Plot the slopes
time_slopes_by_year2 <- time_slopes_by_year %>%
  separate(id, c("X", "Y"), sep = "_") %>%
  mutate(X = as.numeric(X),
         Y = as.numeric(Y))

plot_map_raster +
  geom_raster(data = time_slopes_by_year2, aes(x = X * 1000, y = Y * 1000, fill = estimate)) +
  scale_fill_gradient2(midpoint = 0, name = "Slope (condition factor~year)") +
  geom_sf(size = 0.1)

<img src=“condition_model_cf_files/figure-html/calculate”spatial trend”-1.png” width=“672” style=“display: block; margin: auto;” />


ggsave("figures/supp/condition/spatial_trend_condition.png", width = 6.5, height = 6.5, dpi = 600)

Visualize conditional effects

Depth

# Prepare prediction data frame
nd_dep <- data.frame(depth_sc = seq(min(d$depth_sc), max(d$depth_sc), length.out = 100))

nd_dep <- nd_dep %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model
p_cond_dep <- predict(mfull, newdata = nd_dep, se_fit = TRUE, re_form = NA)

cond_dep <- ggplot(p_cond_dep, aes(depth_sc, exp(est),
                                   ymin = exp(est - 1.96 * est_se),
                                   ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Depth[haul]))

Oxygen

# Prepare prediction data frame
nd_oxy <- data.frame(oxy_sc = seq(min(d$oxy_sc), max(d$oxy_sc), length.out = 100))

nd_oxy <- nd_oxy %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         #oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model
p_cond_oxy <- predict(mfull, newdata = nd_oxy, se_fit = TRUE, re_form = NA)

cond_oxy <- ggplot(p_cond_oxy, aes(oxy_sc, exp(est),
                                   ymin = exp(est - 1.96 * est_se),
                                   ymax = exp(est + 1.96 * est_se)))+
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Oxygen[haul]))

Temperature

# Prepare prediction data frame
nd_tem <- data.frame(temp_sc = seq(min(d$temp_sc), max(d$temp_sc), length.out = 100))

nd_tem <- nd_tem %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         #temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model
p_cond_tem <- predict(mfull, newdata = nd_tem, se_fit = TRUE, re_form = NA)

cond_temp <- ggplot(p_cond_tem, aes(temp_sc, exp(est),
                                    ymin = exp(est - 1.96 * est_se),
                                    ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Temperature[haul]))

Cod

# Prepare prediction data frame
nd_cod <- data.frame(density_cod_sc = seq(min(d$density_cod_sc), max(d$density_cod_sc), length.out = 100))

nd_cod <- nd_cod %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         #density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model
p_cond_cod <- predict(mfull, newdata = nd_cod, se_fit = TRUE, re_form = NA)

cond_cod <- ggplot(p_cond_cod, aes(density_cod_sc, exp(est),
                                   ymin = exp(est - 1.96 * est_se),
                                   ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Cod[haul]))

Sprat

# Prepare prediction data frame
nd_spr <- data.frame(biomass_spr_sd_sc = seq(min(d$biomass_spr_sd_sc), max(d$biomass_spr_sd_sc), length.out = 100))

nd_spr <- nd_spr %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         #biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model
p_cond_spr <- predict(mfull, newdata = nd_spr, se_fit = TRUE, re_form = NA)

cond_spr <- ggplot(p_cond_spr, aes(biomass_spr_sd_sc, exp(est),
                                   ymin = exp(est - 1.96 * est_se),
                                   ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  #coord_cartesian(ylim = c(-4.65, -4.4)) + 
  theme(aspect.ratio = 1) + 
  labs(x = expression(Sprat[sub]))

Saduria

# Prepare prediction data frame
nd_sad <- data.frame(density_saduria_rec_sc = seq(min(d$density_saduria_rec_sc),
                                                  max(d$density_saduria_rec_sc), length.out = 100))

nd_sad <- nd_sad %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         #density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model
p_cond_sad <- predict(mfull, newdata = nd_sad, se_fit = TRUE, re_form = NA)

cond_sad <- ggplot(p_cond_sad, aes(density_saduria_rec_sc, exp(est),
                                   ymin = exp(est - 1.96 * est_se),
                                   ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  #coord_cartesian(ylim = c(-4.75, -4.4)) + 
  theme(aspect.ratio = 1) + 
  labs(x = expression(Saduria~entomon[rec]))

Plot together

# cond_dep + cond_oxy + cond_temp + cond_cod + cond_spr + cond_sad
#   plot_annotation(tag_levels = 'A') 
# 
# ggsave("figures/supp/condition/conditional_effects.png", width = 6.5, height = 6.5, dpi = 600)

p_cond_dep2 <- p_cond_dep %>%
  mutate(variable = "Depth (haul)") %>% 
  dplyr::select(est, est_se, depth_sc, variable) %>% 
  rename(value = depth_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_cond_oxy2 <- p_cond_oxy %>%
  mutate(variable = "Oxygen (haul)") %>% 
  dplyr::select(est, est_se, oxy_sc, variable) %>% 
  rename(value = oxy_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_cond_tem2 <- p_cond_tem %>%
  mutate(variable = "Temperature (haul)") %>% 
  dplyr::select(est, est_se, temp_sc, variable) %>% 
  rename(value = temp_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_cond_cod2 <- p_cond_cod %>%
  mutate(variable = "Cod (haul)") %>% 
  dplyr::select(est, est_se, density_cod_sc, variable) %>% 
  rename(value = density_cod_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_cond_spr2 <- p_cond_spr %>%
  mutate(variable = "Sprat (sub)") %>% 
  dplyr::select(est, est_se, biomass_spr_sd_sc, variable) %>% 
  rename(value = biomass_spr_sd_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_cond_sad2 <- p_cond_sad %>%
  mutate(variable = "S. entomon (rec)") %>%
  dplyr::select(est, est_se, density_saduria_rec_sc, variable) %>% 
  rename(value = density_saduria_rec_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

conds <- bind_rows(p_cond_dep2, p_cond_oxy2, p_cond_tem2, 
                   p_cond_cod2, p_cond_spr2, p_cond_sad2)

ggplot(conds, aes(value, exp(est), ymin = exp(est - 1.96 * est_se),
                  ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  facet_wrap(~variable, scales = "free_x") +
  labs(y = "Le Cren's condition factor",
       x = "Scaled value") +
  theme_plot()


ggsave("figures/supp/condition/conditional_effects.png", width = 6.5, height = 6.5, dpi = 600)